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Abstract 

In  this  paper  we  develop  and  analyze  real-time  and  accurate  filters  for  nonlinear 
filtering  problems  based  on  the  Gaussian  distributions.  We  present  the  systematic 
formulation  of  Gaussian  filters  and  develop  efficient  and  accurate  numerical  integration 
of  the  optimal  filter.  We  also  discuss  the  mixed  Gaussian  filters  in  which  the  conditional 
probability  density  is  approximated  by  the  sum  of  Gaussian  distributions.  A  new 
update  rule  of  weights  for  Gaussian  sum  filters  is  proposed.  Our  numerical  testings 
demonstrate  that  new  filters  significantly  improve  the  extended  Kalman  filter  with  no 
additional  cost  and  the  new  Gaussian  sum  filter  has  a  nearly  optimal  performance. 


1  Introduction 

The  nonlinear  filtering  problem  consists  of  estimating  the  state  of  a  nonlinear  stochastic 
system  from  noisy  observation  data.  The  problem  has  been  the  subject  of  considerable 
research  interest  during  the  past  several  decades  because  it  has  many  significant  applications 
in  science  and  engineering  such  as  navigational  and  guidance  systems,  radar  tracking,  sonar 
ranging  and  satellite  and  airplane  orbit  determination  [11,  14,  15].  As  is  well-known,  the 
most  widely  used  filter  is  the  extended  Kalman  filter  for  nonlinear  filtering  problems.  It  is 
derived  from  the  Kalman  filter  based  on  the  successive  linearization  of  the  signal  process  and 
the  observation  map.  The  extended  Kalman  filter  has  been  successfully  applied  to  numerous 
nonlinear  filtering  problems.  If  nonlinearities  are  significant,  however,  its  performance  can 
be  substantially  improved.  Such  efforts  have  also  been  reported  in  [1,  5,  7,  12,  13,  20].  In 
this  paper  our  objective  is  to  develop  and  analyze  real-time  and  accurate  filters  for  nonlinear 
filtering  algorithms  based  on  Gaussian  distributions.  We  present  the  systematic  formulation 
of  Gaussian  filters  and  mixed  Gaussian  filters. 

We  first  develop  the  Gaussian  filter.  The  proposed  filter  is  based  on  the  steps:  (1)  we  as¬ 
sume  the  conditional  probability  density  to  be  a  Gaussian  distribution  (i.e.,  assumed  density) 
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and  next,  (2)  we  obtain  the  Gaussian  filter  by  equating  the  Bayesian  formula  with  respect 
to  the  first  moment  (mean)  and  the  second  moment  (covariance).  Our  approach  is  based 
on  the  efficient  numerical  integration  of  the  Bayesian  formula  for  optimal  recursive  filtering. 
The  direct  evaluation  of  the  Jacobian  matrix  associated  with  the  extended  Kalman  filter  is 
avoided,  which  is  similar  to  the  one  recently  reported  in  [12].  Secondly,  we  use  Gaussian  sum 
filters  for  the  development  of  nearly  optimal  filters.  The  Gaussian  sum  filter  has  been  studied 
in  [1,  20].  However,  we  adapt  our  Gaussian  filter  for  the  update  of  Gaussian  distributions. 
We  also  suggest  some  new  update  rules  of  weights  for  Gaussian  sum  filters.  Through  our 
experimental  study  we  found  that  the  filters  developed  in  the  paper  perform  better  than  or 
as  good  as  the  filter  of  Julier-Uhlmann  [12].  They  have  a  significant  improvement  over  the 
extended  Kalman  filter  with  no  additional  cost. 

An  outline  of  the  paper  is  as  follows.  In  Section  2  we  develop  the  Gaussian  filter  based 
on  a  single  Gaussian  distribution.  In  Section  3  we  discuss  the  efficient  numerical  integration 
of  the  Gaussian  filter  based  on  quadrature  rules,  and  introduce  the  Gauss- Hermite  filter 
(GHF)  and  the  central  difference  filter  (CDF).  In  Section  4  we  formulate  the  Gauss- Hermite 
filter  and  the  central  difference  filter  as  the  filter  algorithms.  In  Section  5  we  introduce  the 
mixed  Gaussian  filter  and  the  new  update  rules  of  weights.  In  Section  6  we  discuss  the 
relation  of  the  Gaussian  filter  for  the  discrete  time  system  to  the  continuous  time  optimal 
filter  governed  by  the  Zakai  equation.  In  Section  7  we  analyze  the  stability  and  performance 
bound  of  the  Gaussian  filters  and  the  mixed  Gaussian  filters  developed  in  the  paper.  In 
Section  8  we  report  our  numerical  findings  and  comparison  studies.  Also,  we  demonstrate  a 
nearly  optimal  performance  of  the  mixed  Gaussian  filter.  Finally,  we  conclude  our  results  in 
Section  9. 


2  Gaussian  Filters 

We  discuss  the  nonlinear  filtering  problem  for  the  discrete-time  signal  system  for  i?”-valued 
process  x(k): 

(2.1)  x(k)  =  f(x(k  —  1))  +  w(k) 
and  the  observation  process  y(k)  £  Rp  is  given  by 

(2.2)  y(k)  =  h(x(k))  +  v(k) 

where  w(k)  and  r  (  k)  are  white  noises  with  covariances  Q  and  R,  respectively.  We  assume  that 
the  initial  condition  ,7;(0)  and  w(k),  v(k)  are  independent  random  variables.  The  optimal 
nonlinear  filtering  problem  is  to  find  the  conditional  expectation  E(x(k)\Yk)  of  the  process 
x(k )  given  the  observation  data,  V).  =  { y ( j ) ,  I  <  j  <  k}.  The  probability  density  function 
pk\k  of  the  conditional  expectation  E[x(k)  |W]  is  given  by  Bayes’  formula 

(2.3)  ((2,)4,q),ft  dt 
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and 


(2.4) 


where  pk\k-i  is  the  one-step  prediction  and  is  the  probability  density  function  of  x(k)  condi¬ 
tioned  on  Yk_i.  That  is,  the  recursive  filter  (2.3)-(2.4)  consists  of  the  prediction  step  (2.3) 
and  the  correction  step  (2.4).  We  consider  a  Gaussian  approximation  of  this  recursive  for¬ 
mula  and  develop  Gaussian  filters.  We  assume  that  /y-_  i  p-i  is  a  single  Gaussian  distribution 
with  mean  xk_1\k_1  and  covariance  Pk_1\k_1.  Then  we  construct  the  Gaussian  approximation 
of  pk\k  in  the  following  two  steps.  First,  we  consider  the  predictor  step.  We  approximate 
Pk\k-i  by  the  Gaussian  distribution  that  has  the  same  mean  and  covariance  as  pk\k-i-  By 
Fubini’s  theorem  the  mean  and  covariance  of  pk\k- 1  are  given  by 


iM  =  jm  -n>k\k  ,(•'•)  dx  =  JRn 


((2-7T  )ndet  Q)1/2 


rl  k*mYQ-H*-mxdx 


Pk—l\k  —  l  ( f )  dft 


and 


J ^  f  (t'jPk— i|fe— l (^)  d,t 

Ek\k- i[xJ]=  [  xxtpkV^.1(x)dx 

JRn 

=  /R„  (/R„  -mdsow  *,)  it 


Q+f  f{t)f{tYPk-i\k-i{t)dt 

J  Rn 


Remark  1.1:  To  derive  the  mean  and  covariance  of  pk\k-h  we  also  can  use  (2.1)  and  the 
independence  of  w(k)  and  :r(k  —  1).  For  example, 

Ek\k-i\x]  =  E[f(x(k  -  1))  +  w(&)|y*_;i]  =  E[f(x(k  -  1)) |T*_!] 

which  is  precisely  the  same  as  the  above,  though  differently  expressed. 

Thus,  if  Ph~x\k-i  ig  a  Gaussian  with  mean  xk_1\k_1  and  covariance  Pk_1\k_1 ,  then  the 
Gaussian  approximation  of  pk\k-\  has  mean  xk\k_1  and  covariance  Pk\k~\  dehned  by 


(2.5) 
and 

(2.6) 


k\k  —  l 


f  flf ) _ 1 _  p  2  b  a;^-l|fc-l)tG-l|A-l  h  xk-l\k-l)  Jj 

JnJy-U^rdctPk-Mk-i)1/2  e 


Pk\k-1  —  Q  +  %k\k-l)(f(t)  Xk\k-lY  ((2  ir)™det  Pk_uk_i)112 


e  2  h  */.  i|;  i  ,|;  ,  ('  ;r/  i|f  "■  (jj 

Next,  we  discuss  the  corrector  step. 
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where  c  is  the  normalization  constant  and  we  assume  that  pk\k-\  is  given  by  the  Gaussian 
approximation  defined  by  (2.5)-(2.6).  Define  the  innovation  process 

I(k)  =  y(k )  -  Kk\k  i  [// ( k )J  =  y(k )  -  Bk\k-i[h{x(k)]. 


Here  we  again  approximate  for  the  conditional  distribution  of  (x(/s),  h(x(&))  given  the  ob¬ 
servations  up  to  time  A’  —  1  by  a  Gaussian.  That  is,  we  approximate  Ek\k_1[h(x(k)]  by  its 
Gaussian  approximation  z.  This  means  that  the  probability  density  function  of  z  is  given 
by  the  Gaussian  distribution  with  mean  z  and  covariance  Pzz  defined  by 


(2.7) 
and 

(2.8) 


z  = 


f  h(f) _ 1 _ _  p  2  h  xk\k-l)tPk\k_1{t  xk\k-l)  n 


p  = 

zz 


= Jjm  -  mm  - 


Finally,  we  construct  the  Gaussian  approximation  of  pk\k  with  mean  xk\k  and  covariance 
defined  by 


(2.9)  xk\k  =  xk\k_g  +  Lk(y(k )  -  z) 

and 

(2-10)  Pk\k  =  Pk\k-i  ~  LkPl 

where  the  Kalman-hlter  gain  Lk  is  dehned  by 

(2.11)  Lk  =  PXZ(R  +  Pzz)  1 


and  the  covariance  Pxz  is  dehned  by 


(2.12)  Pxz  =  f  (I.  xk\k  ,)(//(/)  z)' 

JR" 


t(,27r) "det  Pk\k_i ) 1  /2 


3  2  h  i*7* |A  — 1 ) I*— 1  (*  'f:*M  1  (It 


In  order  to  implement  the  Gaussian  hlter  we  must  develop  the  approximation  methods 
to  evaluate  integrals  (2.5)-(2.12).  In  Section  3  we  discuss  the  approximation  methods  for 
the  integration  of  the  form 


L  F{t) 


((2  G 


_i _  e“  jf*-1 

ld.et  E)1/2  e 


■Y^-pt-Pdt 


where  ^(A)  is  a  given  function.  Especially  our  discussions  include  the  Gauss-Hermite  quadra¬ 
ture  rule  and  finite  difference  approximation.  Note  that  a  different  approximation  to  this 
integration  results  in  a  variant  of  the  Gaussian  hlter.  For  example,  if  we  approximate  the 
integration  by  linearizing  F  at  x.  then  we  obtain  the  extended  Kalman  hlter. 
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In  order  to  improve  the  performance  of  Gaussian  filter  we  also  discuss  the  mixed  Gaussian 
filter  in  Section  5.  We  approximate  pk~ i|*_i  by  the  linear  combination  of  multiple  Gaussian 
distributions,  i.e. , 


Pk  1 1/;  ](•'•)  =  J2a 

i= 1 


(0 

k- 1 


-L(x-xU-  )-!  (x-P’)  ) 


□  (*.) 


{{2v)ndetPW-  p/2 


(*) 


Then,  each  Gaussian  distribution  is  updated  separately  by  the  Gaussian  filter  (2.5)-(2.12). 
In  the  corrector  step  we  update  the  weights  ak\  For  example,  we  can  determine  the  weights 
Gy-  by  the  T2-projection,  i.e.,  ak\  1  <  i  <  ?u,  minimize 


/ 


a 


(0 


i= 1 


((2rJ"«.  •  /'J  j/  J1  Z* 


to  \ 

k\k> 


dx 


over  (R+)m,  where  pk\k  is  dehned  as  in  the  corrector  step  (2.4),  assuming 


Pk\k- 1(*)  = 


a 


(0 

k- 1 


p(0 


(») 


3  Quadrature  Rules 


f  pa) _ i _  e-§b-<?Ti  ‘0  a 

JRn  r  GJ  ((2'7r)T1(iet  S)1/2  e 


In  this  section  we  discuss  the  approximation  methods  for  the  integral  of  the  form 
(3.1)  /  = 

If  we  assume  S  =  StS  and  change  the  coordinate  of  integration  by  /,  =  Sts  +  x.  then 

<3-2)  1  =  L  P{s)Wr* 


-bP  ds 


e  2 


with  F(s)  =  F(Sts-\-x).  We  apply  the  Gauss-Hermite  quadrature  rule.  The  Gauss- Hermite 
quadrature  rule  is  given  by 


f  OO  2  771  ^ 

/  //W  e_a'  dx  =  X! 

•y  <x  v  ;  i=i 


where  the  equality  holds  for  all  polynomials  of  degree  up  to  2 m  —  1  and  the  quadrature 
points  xt  and  the  weights  are  determined  (e.g.,  see  [9])  as  follows.  Let  ./  be  the  symmetric 
tri-diagonal  matrix  with  zero  diagonals  and  +1  =  yji /2,  1  <  i  <  m  —  1.  Then  {ay}  are  the 
eigenvalues  of  J  and  uy  equal  to  |(r,)i  |2  where  (r.()i  is  the  hrst  element  of  the  i-th  normalized 
eigenvector  of  J .  Thus,  I  is  approximated  by 


(3.3) 


G  =  V  .  .  .  Y  F(qn,qu, ...  ,q, 
1 1  —  1  —1 


,  w. 
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where  qt  =  ^L,  1  <  i  <  m  and  Im  is  exact  for  all  polynomials  of  the  form  snsn  . .  .  stn 
with  1  <  ik  <  2m  —  1.  In  order  to  evaluate  lm  we  need  nfi  point  function  evaluations.  For 
example  m  =  3  we  have 

q1  =  —a/3,  <Z2  =  0,  §3  =  /3  and  uq  =  rc3  =  w2  = 

o  3 

and  J3  requires  9-point  function  evaluations  for  n  =  2. 

In  Julier-Uhlmann  [12]  the  standard  normal  distribution  is  approximated  by  the  discrete 
distribution  as 


(  y/n  +  k  ek,  1  <  k  <n 


Zk  —  <  —Zb_ 


k—n 


n  +  1  <  k  <  2n  and  P(Z k)  = 
k  —  2n  +  1, 


2k  A’  =  2rc  +  1 

1  <  A’  <  2n 


2(tjH-k) 

1 


2(b  +  k)  a  — 


where  k  >  0  is  a  constant  and  ek  is  A’-th  unit  vector  in  R" .  Here,  this  discrete  distribution 
has  the  same  hrst,  second  and  higher  odd  moments  as  the  standard  normal  distribution. 
Julier-Uhlmann  developed  a  Gaussian  filter  based  on  the  following  quadrature  rule  for  (3.1) 


2n+l 

(3.4)  Iju  =  £  F{Zk)P{Zk). 

k=  1 

The  JU-rule  requires  (2 n  +  l)-point  function  evaluation  and  is  exact  for  all  quadratic  poly¬ 
nomials.  If  we  set  k  =  2  and  n  =  1,  then  the  JU-rule  coincides  with  the  Gauss-Hermite  rule 

G- 

Finally  we  consider  the  polynomial  interpolation  methods.  We  approximate  F(t)  by  the 
quadratic  function  P2  that  satisfies 


P{si )  —  P2{si), 


1  <  *  < 


( n  +  1 )(//  +  2) 

2 


at  the  points  {s4}  in  Rn  consisting  of 

0,  h  e* ,  —h  e4,  1  <  i  <  n  and  he4  +  htj.  I  <  i  <  j  <  n, 
where  h  >  0  is  the  stepsize.  That  is,  P2  is  given  by 

(3.5)  -P2(s)  =  P(0)  +  5Z"=1  a*  Sj  +  — 

where  5,  is  the  l-th  coordinate  of  a  point  s  £  if” ,  and  the  symmetric  matrix  PI  =  f  II  j.j  )„Xn 
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are  defined  by 


a,-  = 


F{het)  -  F(-het) 

2/?  ’  “  “ 


n 


II;  ;  = 


H ;  ;  = 


_  —  2P(0)  +  F(  —  hei) 


h 2 


1  <  *  < 


n 


i\h(!j  +  // r; j  j  —  /•  (  hc.j)  —  F(  —  hej )  +  F( 0) 

/P  ’  1 


F  i  <  j  <  n. 


Thus,  we  approximate  (3.1)  by 
(3.6) 

and  the  integral 


h  = 


=  F(0)  +  £?=1 


by 

(3.7) 


Jg  = 


!Rn  P (2-)" 

J  =  d, 

i/hwh2’  WiSF'4"fjs 
=  /J1’/,)2’  +  EE,  (afM21  +  1  #1,14?)  +  E¥)  414? 


where  III  *  is  the  quadratic  approximation  of  /q  for  *  =  1,2. 
Next,  we  consider  the  approximation  Pq  of  F,  dehned  by 


(3.8) 


Pi(»)  =  Ho)  +  e"=i  £lM_ZLiE)  s,  +  E-=1 1  H,t,s] 


which  only  bases  on  the  values  F(dohei)  =  Pi(±he;),  1  <  i  <  n  and  P(0)  =  Pi(0),  and  uses 
the  diagonal  second  order  correction  E”=1  f  Has?  of  the  central  difference  approximation  of 
F  in  (3.5).  Then  the  integrals  /  and  J  are  approximated  by 


(3.9) 

and 

(3.10) 

J,  = 


= 


=  /  Pi  M  w  *  =  m  +  EE,  5  h„ 


Pi(/iei)  -  Pi(-hej)  F2(hei)  -  P2{  hci) 


En 

i= 1 


2  h 


+  4  (o)fi(o)  +  EE,  1 4144 


2  h 

If  we  remove  the  second  order  correction  term  in  (3.9)  (3. 10) ,  then  this  coincides  with  the 
extended  Kalman  filter  with  the  central  difference  approximation  of  the  Jacobian  of  F. 
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4  Filter  Algorithms 


In  this  section  we  describe  the  filter  algorithms  that  are  based  on  the  Gaussian  filter  and 
quadrature  methods  described  in  Section  3.  We  first  consider  the  Gauss- Hermite  filter  based 
on  the  Gaussian-Hermite  quadrature  rule.  Let  ( q, .  ?/■, )  be  the  quadrature  rule  for 


(4.1) 


N 

e~\x\  dx  ~  ^2  wiF{qi ). 
i=% 


Let  .To  and  P0  be  the  starting  values  for  the  mean  and  covariance  of  the  random  variable 
;r(0)  and  set  P0 10  =  Po  and  .r0|o  =  .f0-  In  Gaussian-Hermite  hlter  we  apply  the  following 
predictor  and  corrector  steps  recursively. 

Gauss-Hermite  Filter 

Predictor  Step  Compute  the  factorization  Ph~i\h-i  —  StS  and  set  aq  —  Stc[i  +  xk_1\k_1. 
Update  =  N(xk\k_1,Pk |*_i)  by 

Xk\k-1  =  Ef=l  f(Xi)Wi 

Pk\k—  1  Q  Ei=i  (/(p)  $k\k— 1 ) ( /(-G)  k\k  I  ) 

Corrector  Step  Compute  the  factorization  Pk\k-i  =  TtT  and  set  aq  —  Tf  qi~\-  xk\k_i.  Update 
Pk\k  Ff  (x k\k  i  Pk\k')  t*y 

%k\k  ^k\k  —  l  U  Lk(y(k)  Zf j) 


Pk\k  Pk\k  —  1  PfcPxz 

where 

zk  =  Th= i  h(xi)wi 

Pxz  E*=1  (*G  ^k\k—  1 ) zk)  i ■ 

Pzz  =  Ejll  (MX«)  -  -  zkfwi. 


Lk  —  Pxz(P  +  Pzz }  '  • 

We  use  the  Choresky  decomposition  for  finding  S  and  T  [6]  in  our  calculations.  One  of 
the  advantages  of  the  quadrature  based  hlter  is  that  we  do  not  require  to  have  the  derivative 
of  f  and  h. 

Secondly,  we  present  the  hlter  algorithm  based  on  the  central  difference  approximation 
with  second  order  diagonal  correction.  We  summarize  the  hlter  algorithm  of  the  central 
difference  approximation  with  second  order  diagonal  correction  as  follows. 

Central  Difference  Filter 
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Predictor  Step  Compute  the  factorization  Pk-i\k-i  —  StS  and  set  f(s )  =  f[Sts-\-xk_ i|fe-a) 
approximated  by  /( 0)  +  £"=i  atst  +  \ s^s.  Update  =  N(xk\k_1,  Pk\k-t^  by 

xk\k- 1  =  /(O)  +  Y!i= 1  2 

n-p-1  =  Q  +  E?=1  +  E?=1  ^ 

where  H \  ,■  are  the  central  difference  approximations  of  'J  . 

csi 

Corrector  Step  Compute  the  factorization  =  TfT  and  set  h(s)  =  h(TU  + 

approximated  by  h(0)  +  E”=i  G'5?  +  \stGs.  Update  pk\k  =  ■\'(xk\k.  Pk\k)  by 

—  xk\k—i  T  Lk(y (h)  zk^j 

Pk\k  Pk\k  —  1  LkPxz 

where 

**  =  m  +  i  Gw 

Cr  =  '/''(/'I  -  •  •  •  :  Kf- 

P„  =  E”=1M'  +  E"=1iG,„G',, 

Lk  —  PUR  +  P?.r  )  '  • 

where  Gj,j  are  the  central  difference  approximations  of  0  . 

In  the  algorithm  we  avoid  to  calculate  the  derivatives  of  /  and  h.  Instead,  we  use  the 
central  difference. 

The  followings  are  the  operation  counts  in  terms  of  function  evaluations. 

If  we  use  the  quadrature  rule  based  on  the  g-point  Gauss- Hermite  rule  then  the  algorithm 
requires  (n  +  p)qn  function  evaluations. 

If  we  use  the  Julier-Uhlmann  discrete  Gaussian  rule  then  the  algorithm  recpiires  (n+p)(2n+l) 
function  evaluations. 

If  we  use  the  central  difference  algorithm,  the  (n+p)(2n  +  l)  function  evaluation  are  required. 

As  we  will  see,  our  numerical  testings  indicate  that  the  algorithm  based  on  the  Gauss- 
Hermite  rule  performs  better  than  the  others.  If  we  use  the  3-point  rule,  then  we  require 
729(n  +  p)  function  evaluations  for  the  case  n  =  6,  which  can  be  done  in  real-time. 
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5  Mixed  Gaussian  Filter 


In  this  section  we  discuss  the  mixed  Gaussian  filter.  We  approximate  the  conditional  prob¬ 
ability  density  pk~\ p_!  by  the  linear  combination  of  multiple  Gaussian  distributions,  i.e., 


(5.1)  pk  ,|A.  ,(.r)  = 


JO 


2  '  fc-llA-l^  '  A-llA-1'  N?  ^A-llA-l' 


□(«) 


a"_1  (WS^IA-J172 


03 


Here  we  apply  the  Gaussian  filter  (2.5)-(2.12)  to  each  Gaussian  distribution  fV(x^1|A_1,  JP|^1|i._1) 

and  obtain  the  update  N(x^k,Pf^k).  Each  update  is  independent  from  the  others  and  can 
be  performed  in  a  parallel  manner. 

Next,  we  update  the  weights  aj^  for  the  new  update  pk\k{x )  at  the  end  of  corrector  step. 

We  discuss  here  three  update  formulae  in  what  follows. 

First,  by  equating  the  first  moment  of  each  Gaussian  distribution  we  obtain 


*:<»• 


.(«)  ■) 

A|A;  ,/./• 


- J 


=  a 


.(0  / 

‘-1  h 


(y-Kx))tR  )*(**!*-! )  h^-^AiA-i))  dx 


Rn  (2'K)n{detRdetPk^k_1) 


<'.)  .1/2 


Here  we  approximate  the  right  hand  side  by  the  Gaussian  distribution  as  in  (2.5)  (2.6)  and 
obtain 


(5.2) 


JO  _  JO 


_ 1  -ife-DJR+Fj- J-i) 

/'•  I  ((2ir)ndet  (R+P^))1/2 


where  £,  Pzz  are  dehned  by  (2.7)— (2.8),  which  is  the  update  formula  discussed  in  [1], 


Next,  we  apply  the  collocation  condition  at  x 


(0 

k\k  * 


a 


(0 


({'2r)rvi,:<  Pj,  )l/2 


_  JO 


=  a 


fc— l ' 


=- hi  M^i;  -4;Li )  JPA|A-1 )  1  Gaia  :r;  i;  0) 


((27r)n<ie#W_1)1/2  ” 

to  obtain  the  update 

(5  3)  a(0  _  (0  Gef4iA>1/2  -|(h-M^A))tfi~J^-fe(4|A))+(4|A-4|L1)t(FAlA-1)~1(4'l~a;^-i)) 

1  J  k  k~x  (^4;Li)1/2 

Finally  we  discuss  the  simultaneous  update  of  the  weights.  We  determine  the  weights 
ak'  by  the  P2-projection,  i.e.,  1  <  i  <  m  minimize 


(5.4) 


/ 

JR" 


i= 1 


JO _ 1 _ 

k  {(2^detP^y/l 


Pk\k{x )  /-M..,.  e 


U/p(0  )-l  (x— J’  )  \ 
2  .  •  ;  |;  •  •  ;  i;  • 


a  |  a  ' 


dx 
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over  (R+)m:  where  p f.\k  is  defined  as  in  the  corrector  step  (2.4)  with 


Pk\k-i(x)  =  E 


a 


(0 
k- 1 


((2-)"^iiLi)i/2 


-i) 


In  order  to  perform  the  minimization  (5.4)  we  need  to  evaluate  the  integral  of  the  form 

JRne  ((2^detXyn  e  aX- 

and  it  is  relatively  expensive.  Hence  we  propose  the  minimization  of  the  sum  of  collocation 
distances: 


(5.5) 


E 

8  =  1 


w  t  _  v-  Ji) _ 1 _ 


P*l*(a:*|*)  -  E 


a 


j=i 


k  ({2^detP[\l)U2 


over  a  6  satisfying  o  >  a0  >  0.  A  positive  constant  o0  is  chosen  so  that  the  likelihood 
of  each  Gaussian  distribution  is  nonzero  (e.g.,  a0  =  0.001(1, .  .  .  ,  1)*  £  Rm).  Problem  (5.5)  is 
formulated  as  the  quadratic  programming 


(6.6) 


min  —  atAtAa 


subject  to  a  >  a o, 


where  8  >  0  is  chosen  so  that  the  singularity  of  the  matrix  A 1 A  is  avoided  and  the  matrices 
(A,  b)  are  defined  by 


Aii  — 


((2*  rdetP^yn 


40 

k\k 


J.f)  ’ 

k\  k> 


a(0  )—  1 1 

1  k\k> 


4  0  _r(j)  ) 

'  k\k  k\k' 


b*  =  E 


(2^(detRdetP^[_1P^ 


))ti?  1#-/l(4'|i))+(3 


.(•)_„(«)  it/pO.)  6-1/. 
fc|A  ^ 7? | /?  —  1  /  Wfclfc-I.)  V 


r(J) 

^  fc|A  - 


-a- 


Thus,  we  solve  (5.6)  to  obtain  the  weights  ajp  at  each  corrector  step  by  using  the  existing 
numerical  optimization  method  (e.g.,  see  [6]). 

The  theoretical  foundation  of  the  Gaussian  sum  approximation  as  above  is  that  any 
probability  density  function  can  be  approximated  as  closely  as  desired  by  a  Gaussian  sum. 
More  precisionly,  we  state  the  following  error  estimate,  the  proof  of  which  is  found  in  [17]. 

Theorem  5.1:  Let  A4  be  a  non- negative  integer  and  M  =  2 A4  +  2.  For  any  e  >  0  there 
must  exist  D  >  0  and  a  mask 


l  j  =  (ji,  ■  ■  ■  jnf  g  zn,  \j\  =  jr  js  <  x| 

such  that  for  any  density  function  p[x  )  e  C^{Rn)  n  W*{Rn)  and  all  h  >  0  the  estimate 
(5.7)  II p(x)  -  Ph(x)  II  <  a-/'A’|pkv.xv'(/."']  +  (Wl>Ww£  V>G 
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holds,  where  pn(x )  is  a  linear  combination  of  Gaussian  distributions  given  by 

1 


(6.8) 

with 


Ph{x)  = 


y^TT )nDh2  jfcft 


_  \x-hk\2 

2^  !>k  <  2Dh2 


Pk 


=  E  n,  -  hi) 


2  OLA 


\j\<M 


and  cy  is  a  constant  independent  ofp(s).  Here  |  •  and 

\p\w£(R”)  =  E  ll5‘?Pllc(i?") 

\3\=M 


are  defined  by 


and 


Af-1 


E 

lil=o 


\&p\ 

j! 


respectively. 

Roughly  speaking,  the  estimate  (5.7)  shows  that  any  probability  density  function  p(x)  £ 
CAf(Rn)  fl  W^(Rn)  can  be  approximated  by  a  sum  of  Gaussian  distributions  each  of  whose 
components  is  given  by 


J  \x—hk\* 

P  2  D  h2 

(27t)™  Dh 2 

with  order  0{h^)  for  h  — *  0. 


6  Relation  to  Continuous-Time  Filter 

In  this  section  let  us  discuss  the  nonlinear  filtering  problem  for  the  continuous-time  signal 
process  x(t)  in  Rn  generated  by 

(6.1)  dx(t)  =  /(:/•(/.))  dt  +  Q1/2dB(t). 

where  B(t)  is  the  standard  Brownian  motion,  i.e. ,  x(t )  is  the  diffusion  process  [3,  11,  19]. 
So,  equation  (6.1)  holds  in  the  sense  of  Ito  and  x(t)  satisfies 

x(t)  =  s(s)  +  f  f(x(T))dT  +  Q1/2(B(t)  -  B(s)) 

J  S 

for  all  0  <  s  <  t. 

Now  we  consider  the  continuous  observation  process  z{t)  £  Rp 

(6.2)  z(t)  —  j  h(x(s))  ds  +  R1/2V(t) 


12 


where  V  (t)  is  the  standard  Brownian  process  that  is  independent  of  B f  /  j  and  x  (0). 

Or,  we  consider  the  discrete  observation  process  y(k)  £  Rp 

(6.3)  y(k)  =  h(x(kAt))  +  v(k). 

where  At  >  0  is  the  stepsize  and  v(k)  is  white  noise  with  covariance  R.  We  assume  that  the 
initial  condition  x(0)  and  B(t ),  v[k)  are  independent. 

Then  the  conditional  probability  density  p(t)  =  p(t:x )  £  L2(Rn )  satisfies  the  Zakai 
equation 

(6.4)  dp(t)  =  Ap(t)dt  +  R~1hpdz[t ) 

where  A  is  the  Fokker-Planck  operator  [3,  19]. 

As  shown  in  [8,  10,  15],  the  discrete-time  filter  (2.3)-(2.4)  applied  to  the  time-discretized 
signal  system  of  (6.1)-(6.2): 

x(k)  =  F(x(k  -  1))  +  Q1/2(B{kAt )  -  B((k  -  1)A t)) 

(6.5) 

y(k)  =  H(x(k ))  +  R1/2(V(kAt)  -  V((k  -  1)A t)) 

provides  an  approximation  method  for  the  continuous-time  optimal  filter  to  (6.1)-(6.2).  Note 
that  F(s )  denotes  the  finite  difference  approximation 

x  ~  F{s)  =  s  +  At  f(s) 

for  the  ordinary  differential  equation: 

=  f  {■<'').  ;r(0)  =  5,  .S  £  Rn. 

Next,  we  apply  the  discrete-time  filter  in  the  following  manner.  We  discuss  the  continuous¬ 
time  process  (6.1)  with  the  discrete-time  observation  (6.3).  We  approximate  the  process  x{t) 
by  the  discrete-time  process  xJk  by 

(6.6)  x\  -  $(h,  x3-1)  +  Q1/2wJk,  1  <j<  M 

with 

wJk  =  B  (k  At  +  j  5)  —  B  [k  At  +  ( j  —  1 )  6) . 

Here  the  interval  [A’ At,  ( k  +  l)At]  is  subdivided  into  the  M  subintervals 

[kAt  +  (j  —  1)<5,  kAt  +  j8] 

where  1  <  y  <  M  and  8  —  jj-. 

So,  xk  denotes  the  approximation  of  the  process  x{t )  at  t  =  kAt  +  j8  and  $(<5,  .x).-1  j  denotes 
the  discretization  of  the  ordinary  differential  equation 
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We  employ  a  family  of  Runge-Kutta  methods  in  our  calculation,  including  the  Euler  approx¬ 
imation: 


®{s,4  x)  =  4  1  +sf(4  x) 

The  predictor  step  (2.5)-(2.6)  is  applied  successively  M  times  to  obtain  [xk\k_1:  Pk\k_ 1),  he., 


(6.7) 

and 


x 


h\k— 1 


L  m 

JRn 


p3~ 1 
k\k—l 


°k\k 


-i) 


dt 


pj 

rk  |ft-l 


(6.8) 


Q  +  /«»(/«  -  4v,-i)W)  -  4,k-i 


Xe 


-W-xi\k- 


pj— 1 
k\k-l 


h— J_1 


« |A) 


-i) 


dt 


with 


0  _  73  0  _  73 

xk\k-l  ~  xk-l\k-li  rk\k-l  ~  *k-l\k-l 

for  I  <  j  <  M .  Hence,  the  mean  and  covariance  of  pk\k-i  are  given  by 

(6-9)  Xk\k- 1  =  4\k  l:  Pk\k-1  =  Pk\k- 1- 

Then  the  corrector  step  (2.7)-(2.12)  is  applied  as  it  is  to  obtain  (xk\k,  pk\k )■ 


7  Stability  Analysis  of  Filter  Algorithms 

In  this  section  we  analyze  the  stability  and  performance  bound  of  the  Gaussian  filters  (2.5)- 
(2.12)  for  the  signal  process  (6.5)  with  F(s)  =  s  +  A t  f(s),  i.e.,  a  discrete  approximation  of 
the  continuous  signal  system  (6.1)  can  be  given  by 

x(k)  =  x(k  —  1)  +  A t  f(x(k  —  1))  +  \fK tw(k) 

where  w(k)  is  the  white  noise  with  covariance  Q. 

We  assume  that 

(7.1)  \f(xi)  ~  /(®2:),au  -  x2\ <  io\xi  ~  ^2 12,  xh  x2  G  pn ■ 

Then  we  have  the  following  theorem. 

Theorem  7.1:  The  solution  Pk\k  to  (2.5)  (2. 1 2)  has  the  estimate 

(7.2)  trPk\k  <  (1  +uAt)2k  ^rP0|o+  ^  +  J^ty  _  ^  Qj  >  k>0. 


14 


Proof:  From  (2.5)-(2.6)  the  Gaussian  approximation  of  pk\k~ 1  has  the  mean  Xk\k- 1  and 
covariance  Pk\k-i  dehned  by 


(7.3) 


%k\k— i  /  F (s)pk—i^k—i  ds 
JR™ 


and 

(7.4)  iV, 

Note  that 


At  Q  + 


L  M  — l)  X).\k  —  1^  Pk  —  1|A  —  1  (^)  ds. 


F(t )  -  Zfcijfc— !  =  l  (F(t)  -  F(s ))  pk  ||/,.  i(.s)  c/s 


=  /  (/  -  S  +  At(f(t)  -  f(s)))pk_llk_1(s)ds 
JR™ 

Thus  by  the  Cauchy  and  Schwarz  inequality  we  have 


tr 


(m 


^k\k—l 


(m 


%k\k  —  1 


<  /  \t-s  +  At(f(t)-f(s))\2pk_1^1(s)ds. 

JR™ 

Hence  it  follows  from  (7.3)  that 

tr  Pk\k~i  =  At  (tr  Q)  +  tr  jf  ^  jf  ^  \t  -  s  +  At(f(t)  -  f(s))  |2  pJt_i|A:_i(,s)pA,_1|A,_1(/)  ds 


dt 


<  At  (tr  Q)  +trPk_llk_1  +  2AtutrPk_1\k_1  +  Ai2w2hPM|M 
by  assumption  (7.1).  That  is, 

(7.5)  tr  Pk\k_j  <  (l  +  2A t  td  +  At 2  u;2)  tr  +  At  (tr  Q) 

Note  that  for  any  k  the  covariance  Pk\k  is  nonegative  and  satisfies  the  following  property 

Pk\k  —  Pk\k-1- 

Thus  by  (7.5)  we  obtain  the  estimate  (7.2).  □ 

Remark  7.1:  Using  the  similar  arguments  that  lead  to  (7.2),  we  can  prove  that  the  same 
result  holds  for  the  approximation  method  discussed  in  Sections  3-4.  That  is,  the  estimate 
(7.2)  is  still  valid  for  the  Gaussian-Hermite  filter,  provided  that  JZiLi  wi  —  G  wi  >  0. 

Next,  we  discuss  the  stability  of  the  mixed  Gaussian  filter.  The  mixed  Gaussian  filter 
assumes  the  form 

m 

(7-6)  pk\k{x)  =  J2  lAk\k){?) 

.1=1 
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where  the  j-th  Gaussian  distribution  N(x^k,P^)(x)  is  computed  by  the  Gaussian  filters 

described  in  Sections  3-4  in  a  parallel  manner.  The  weights  .  I  <  j  <  rn  are  determined 
at  the  end  of  the  corrector  step  based  on  the  likelihood  of  each  Gaussian  distribution.  Under 
the  assumption  (7.1)  each  Gaussian  filter  in  the  Gaussian  sum  (7.6)  is  stable  in  the  sense  of 

(7-2). 

Next  we  discuss  the  performance  bound.  In  general,  let  the  operators  Vk\k-i{-)  and 
®*|fc-i(')  denote  the  one-step  update  of  conditional  probability  density  functions  that  are 
based  on  the  exact  filter  (2.3)-(2.4)  and  the  mixed  Gaussian  filter  discussed  in  Section  5, 
respectively.  We  estimate  the  difference  between  the  exact  filter  and  the  mixed  Gaussian 
filter: 


(7.7) 


£n  =  ttnk=lVk\k-l  (Po{x))  ~  [  Y  aA $(*) 

U=1 


where  po(x)  is  the  probability  density  function  of  random  variable  x(0)  and 


Po\o  =  Y  aiPo\ 
j= i 


U) 

o 


is  the  mixed  Gaussian  approximation  of  po(x).  Let  |/|  =  fR„  |/|  dx.  Then  we  have  the 
following  result. 

Theorem  7.2:  Assume  that  for  any  probability  density  function  p(x)  and  q(x) 


(7.8) 

Then  the  error  estimate 


(7.9) 


«<7mi  +  D7’-V. 

i=i 


holds,  where  £0  is  the  initial  error  defined  by 


and  the  error  £j-\  is  given  by 


£j- 1  — 


£o  =  Po(x)  -  Y  ajPo\oix} 
3  =  1 


PjU  i(/h  i(-'G  )  -  vp.;i.;  i(/h  i(-T)) 


(7.10) 

Proof:  Note  that 

£n  ~  'Pnln-l  'P’k\k-1  {‘PO  (^)))  —  ^  n\n-l  (Po(®)) 


+  VJ/n|n-l  (n)(=i1Pfc|jl._i(po(^)))  —  ^n\n-l  |  I7fe=l  — 1  (  Y  aj  PO\o{X)  ) 

3  =  1 
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Thus  we  have 

\^n\  ^  7  \^n— 1  |  T  &n—  li 

which  leads  to  (7.9). 

Remark  7.2:  Suppose  po(x)  is  a  mixed  Gaussian  distribution  and  we  set  poio^)  =  po(x). 
Also  assume  that  f(x )  and  h\x)  are  affine  functions.  Then  we  can  show  that  \£n\  =  0, 
i.e.,  the  mixed  Gaussian  filter  is  exact,  since  the  Gaussian  filter  is  exact  for  the  Gaussian 
distribution  in  this  case. 

To  understand  the  validatity  and  usage  of  error  formula  (7.9)  we  adopt  the  following  filter 
(which  is  different  from  the  one  we  implemented),  i.e.,  the  mixed  Gaussian  filter  with  fixed 
means  and  covariances.  For  simplicity  of  our  discussion  we  assume  n  —  1.  Let  us  consider 
the  fixed  Gaussian  elements: 

1  fo-xjf 

=  — ,  e  2<t  1  <  7  <  m 

with  nodal  points  Xj  £  R  and  covariance  u  >  0. 

Then  we  define  the  subspace  Vm  of  the  probability  distributions  on  R  by 

Vm  =  span  ■ 

Dehne  Pm  is  the  projection  of  L2(R )  onto  Vm. 

We  consider  the  mixed  Gaussian  filter  with  hxed  means  and  covariances,  i.e., 

Pk(x)  =  pm  (vj/A| k-iPk-i)  =  \{pk  i ),  Pk-i  e  vm 

where 

m 

Pk{x)  =  e  vrn- 

j= i 

Thus  is  the  hxed  operator  in  this  case.  Hence  £j,  j  >  0  can  be  estimated  based  on 

the  error  estimate  formula  (7.9),  assuming  the  smoothness  of  pj_1[x). 

8  Numerical  Results  and  Comparisons 

In  this  section  we  demonstrate  the  feasibility  of  our  proposed  filter  algorithms,  Gauss- Hermite 
filter  (GHF)  and  central  difference  filter  (CDF),  using  the  three  test  examples  as  in  Examples 
8. 1-8.3.  We  compare  the  performance  of  our  proposed  filters  against  that  of  the  extended 
Kalman  filter  (EKF)  and  the  filter  of  Julier-Uhlmann  (JUF)  [12], 

We  use  the  average  root  mean  square  error  for  our  comparison  of  the  methods.  The 
average  root  mean  square  error  is  defined  by 

(8.1)  GW  = 
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which  is  based  on  N  different  simulation  runs.  Here  subscript  j  denotes  the  j-th  component 
of  the  vector  x(kj  and  its  corresponding  estimate  x(k),  and  the  superscript  i  denotes  the 
2-th  simulation  run. 

Note  that  for  one  dimensional  system  the  filter  of  Julier-Uhlmann  [12]  coincides  with 
the  3-point  Gauss-Hermite  hlter  when  k  —  2.  Also  when  the  stepsize  is  \/3.  the  central 
difference  hlter  has  the  same  estimate  algorithm  as  the  hlter  of  Julier-Uhlmann  except  for 
the  estimation  Pxz.  Hence,  in  this  section  the  stepsize  is  always  chosen  as  \/3  for  our  study. 

We  also  note  that  our  proposed  hlter  algorithms  maintain  the  non-negative  definiteness 
of  the  covariance  update. 

We  carried  out  our  comparison  study  using  various  examples  and  different  starting  states 
a;0|o  and  different  problem  data,  but  we  only  present  the  selected  results  in  what  follows. 

Example  8.1  We  consider  the  one-dimensional  signal  process 

(8.2)  x(k)  =  x(k  —  1)  +  A t  f(x(k  —  1))  +  w(k) 

and  the  observation  process 


(8.3) 


y{k)  =  At  h(x(k ))  +  v(k) 


where 

/'(:r)  =  hx{\  —  x2),  h(x)  =  (x  —  0.05)2 

and 

w{k )  ~  N( 0,  b2At)  v(k)  ~  2V(0,  d2 At) 

It  is  motivated  from  the  continuous  time  systems  (6.1)  and  (6.2)  of  the  form 


(8.4) 


dx(t)  =  f(x)dt  +  bdw(t) 
dy(t)  =  h[x)dt  +  ddv(t) 


as  described  in  Section  6. 

This  example  is  used  as  a  benchmark  because  of  the  following  reasons.  The  deterministic 
equation 


(8.5)  x  =  f(x) 

has  two  stable  equilibria  -1  and  1,  and  one  unstable  equilibrium  0.  So,  we  observed  that  the 
signal  process  X}.  is  distributed  around  one  of  stable  equilibria.  The  observation  function 
h(x)  —  (x  —  0.05)2  is  the  scjuare  of  shifted  distance  from  the  origin  and  distinguishs  the  two 
stable  equilibria  marginally.  The  performance  of  the  hlter  is  possibly  tested  because  of  the 
minor  margin  of  detecting  measurement  function  /t(x)  of  the  signal  process. 

We  use  the  following  problem  data:  the  time  elapse  At  =  0.01;  the  initial  condition 
x0  =  —0.2,  i.e. ,  ,r(0)  =  —0.2;  the  starting  point  h0|0  =  0.8,  P0|o  =  Po  =  2;  the  system 
parameters  b  =  0.5  and  d  —  0.1.  We  consider  the  time  from  0  to  4.  k  is  chosen  as  the  best 
value,  i.e.,  n  =  2  (see  [12]), 
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One  dimensional  system 


Figure  1:  The  Average  of  Root  Mean  Square  Error 


Figure  1  shows  the  average  root  mean  square  errors  committed  by  each  filter  across  a 
simulation  consisting  of  50  runs.  From  Figure  1  we  see  that  the  5-point  Gauss-Hermite  filter 
performs  better  than  others.  The  central  difference  filter  has  similar  performance  like  the 
filter  of  Julier-Uhlmann  [12]  in  this  example.  All  the  three  filters  perform  superior  to  the 
extended  Kalman  filter. 

Example  8.2  In  this  example,  we  consider  the  Lorenz  system 

x(k )  =  x(k  —  1)  +  At  f(x(k  —  1))  +  bw(k),  y(k )  =  At  h{x(k))  +  dv(k ) 

where  x  =  [aq,  x2,  £3]*  is  the  A3- valued  process.  The  if-valued  processes  w(k)  and  v(k)  are 
white  noises  with  the  same  covariance  At.  Here  b  £  A3  is  a  constant  vector  and  d  £  R  is  a 
constant  scalar.  Also 

(8.6)  f(x)  =  [r1(-x1  +  x2),  r2x1  -  x2  -  xxx3,  -r3x 3  +  aqaqf 

The  observation  function  is  chosen  as  the  shifted  distance  from  the  origin  given  by 

h  =  \Jix-x  —  0.5)2  +  x2  +  x\. 

It  is  motivated  from  the  Lorenz  stochastic  differential  system  of  the  form  (8.4)  as  described 
in  Section  6  where  x  £  A3. 

The  three  dimensional  deterministic  equations 

(8.7)  x  =  fix) 
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Lorenz  system  (Component  1 ) 


Figure  2:  The  Average  of  Root  Mean  Square  Error 

is  the  Lorenz  equations  in  which  jq  (  j  =  1,2,3)  are  positive  parameters.  The  three  param¬ 
eters  have  a  great  deal  of  impact  on  the  system.  Here,  we  choose  rq  =  10,  r2  =  28,  and 
r3  =  8/3  which  is  a  mathematically  very  interesting  case,  because  in  the  case  there  are  three 
unstable  equilibria  and  a  strange  attractor  for  the  equations. 

We  use  the  following  problem  data:  the  system  parameters  are  chosen  as  b  =  [0,  0,  0.5f 
and  d  =  0.065;  the  initial  condition  is  x0  =  [— 0.2,  —  0.3,  —  0.5]*,  i.e.,  x(0)  =  [— 0.2,  —  0.3,  —  0.5]*, 
and  At  =  0.01.  We  consider  the  time  from  0  to  4.  The  initial  estimate  is  .r0|0  =  [1.35,  —3,  6 f 
with  covariance  P0|o  =  0.35/. 

In  Figures  2-4  we  show  the  average  of  root  mean  square  errors  for  each  component  of 
system  state  committed  by  the  algorithms  of  the  Gauss-Hermite  quadrature  rule,  the  central 
difference  approximation  and  the  JU-rule  respectively  across  a  simulation  consisting  of  50 
runs.  As  shown  in  Figures  2-4,  the  GHF  and  the  CDF  have  smaller  errors  than  the  JUF  and 
the  EKF.  From  these  figures  we  conclude  that  in  this  example  the  GHF  has  a  substantial 
improved  performance,  and  the  GHF  performs  a  little  better  than  the  CDF. 

Example  8.3  We  further  discuss  the  three-dimensional  continuous  signal  process 

(8.8)  x(t)  =  f(x(t))  +  w(t) 
and  the  discrete  observation  process 

(8.9)  y(t)  =  \J(M2  +  [aq (t)  -  H ]2)  +  r(t) 

where  w(t)  £  R3  and  r f  / )  £  R  are  zero  mean  and  uncorrelated  noises  with  covariances  given 
by  Q(t)  £  R3x3  and  R(t )  6  R,  respectively.  The  function  f(x)  is  given  by 

f(x)  =  [-X2,  -e~7:riX2X3,  of. 
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Average  root  mean  square  error  Average  root  mean  square  error 


Lorenz  system  (Component  2) 


Figure  3:  The  Average  of  Root  Mean  Square  Error 


Lorenz  system  (Component  3) 


Figure  4:  The  Average  of  Root  Mean  Square  Error 


The  above  system  states  x} .  x2  and  x3  represent  altitude  (in  feet),  velocity  (in  feet  per 
second)  and  constant  ballistic  coefficient  (in  per  second),  respectively.  The  detail  physical 
meaning  of  the  system  and  its  parameters  can  be  found  in  [2], 

We  choose  this  example  as  a  benchmark  because  it  contains  significant  nonlinearities  in 
the  signal  and  observation  processes,  and  had  been  discussed  widely  in  the  literature. 

In  the  previous  literature,  a  fourth-order  Runge-Kutta  scheme  with  64  steps  between 
each  observation  is  employed  to  solve  numerical  intergration  of  (8.8)  based  on  the  significant 
signal  process  [2,  12],  That  is,  in  the  predicator  steps  of  EKF  and  JUF,  the  means  are 
calculated  using  the  numerical  scheme.  Then  their  covariance  are  propagated  from  the  jth 
to  (j  +  1 ) t h  step  using 

=  +  ( j  +  1  )<M*  •  )  !’>,-]  fi'tj/  (I.  k  +  ( j  +  1  )<M&  +  j6) 

with 

^{tk  +  (j  +  1)<5,  tk  +  jfi)  —  I  +  SJf  +  2<S’2"^/ 

where  J  j  was  evaluated  at  t +  nS  and  8  =  At/64.  Thus,  P/.+i\k  =  Ptk  +645  \tk- 

In  the  comparison,  we  implement  the  proposed  approximation  schemes  (6.7)-(6.9)  to 
solve  the  hltering  problem  of  (8.8)-(8.9).  In  the  predictor  step  of  GHF,  we  directly  use  the 
Euler  approximation  of  signal  process  (8.8)  based  on  32  or  64  steps  (i.e.,  M  =  32  or  64) 
between  each  observation  (they  are  shortly  denoted  by  GHF32  or  GHF64,  respectively). 
However,  in  NGHF  the  signal  process  system  (8.8)  is  rewritten  as  the  discrete  form 

=  +  $F{xi~1)  +  VSw 

where  x\  is  the  approximation  of  the  process  x{t)  at  t  =  k  At  +  j  6,  and  function  F(s )  = 
(jr (A’i  T  2A’2  T  2k3  T  A^)  with 

A’i  =  f(s,  8),  k2  =  f(s  +  0.5<5 A’!,  S), 


k3  —  f{s  +  0.5hA^2,  ^),  k-4  —  f(s  +  8k3,  8), 


which  is  the  fourth-order  Runge-Kutta  scheme.  Especially,  8  =  AA/64  in  NGHF64. 

We  use  the  following  data:  the  system  parameters  are  chosen  as  7  =  5  X  1 0 — 0 ,  PI  =  10°, 
M  =  10s,  Q(k)  =  0  and  R[k)  =  1 04 ;  the  initial  condition  is  a?(0)  =  [3  X  1 05 ,  2  X  1 0— 4 ,  10-3]4 
and  At  =  0.1.  We  consider  the  time  from  0  to  30.  The  initial  estimate  is  x0|0  =  [3  X  10s,  2  X 
10— 4 ,  3  X  10-5]4  with  covariance 


-Po|o 


106  0  0 

0  4  x  106  0 

0  0  10“4 


We  also  choose  the  optimization  number  k  =  0  in  JUF  for  the  three  dimensional  system. 

In  Figures  5-7  we  show  the  absolute  value  of  average  error  for  each  component  committed 
by  the  algorithms  of  the  Gauss-Hermite  quadrature  rule,  the  JU-rule  and  the  EKF  respec¬ 
tively  across  a  simulation  consisting  of  50  runs.  These  figures  show  that  in  this  example  the 
GHF  has  superior  performance  than  the  JUF  and  the  EKF. 
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Time  (sec) 


Figure  5:  The  Absolute  Value  of  Average  Altitude  Error 


Time  (sec) 


Figure  6:  The  Absolute  Value  of  Average  Velocity  Error 


x  10-3 


Time  (sec) 

Figure  7:  The  Absolute  Value  of  Average  .r3  Krror 

8.1  Mixed  Gaussian  Example 

In  this  section  we  hrst  demonstrate  a  nearly  optimal  performance  of  the  mixed  Gaussian 
method  described  in  Section  5.  We  consider  the  same  example  as  in  Example  8.1  with 
h(x)  —  \x  —  .5|.  We  employ  the  two  Gaussian  distributions  starting  from 

A1)  _  r  pi1)  —  i  aTiri  T(2)  _  _  c  pi2)  —  i 

i(o|o  —  -'T  Wlo  —  -1  dllu  i(o|o  —  -A  ^oio  — 

respectively  and  the  initial  weights  are  q-q1*  =  Oq2'  =  .5.  In  order  to  demonstrate  the 
effectiveness  of  the  proposed  update  (5.5)-(5.6)  we  jumped  the  process  x %  at  k  =  150  (which 
corresponds  to  t  =  1.5  in  the  continuous  time  process  (8.4))  to  x150  =  .5.  It  may  represent  an 
impulse  force  at  t  =  1.5.  We  compare  the  performance  of  the  mixed  Gaussian  filter  against 
the  one  of  the  optimal  Zakai  filter  (6.4)  with  At  =  .01  in  Figures  8-11.  Here  GHF1  (GHF2) 
represents  the  3-point  Gauss- Hermite  filter  which  is  the  hrst  (second)  component  of  the  two 
mixed  Gaussian  filter  with  Weight  1  (Weight  2).  We  approximate  the  Zakai  equation  by 
the  operator-splitting  method  as  described  in  [10].  We  observed  that  the  mixed  Gaussian 
filter  performs  nearly  optimally  up  lo  /  =  1.5.  We  also  note  that  the  Zakai  filter  is  no  longer 
optimal  after  k  =  150  because  of  the  jump.  As  seen  as  Figures  8-9,  the  update  formula 
(5.5)-(5.6)  quickly  captures  the  phase  change  from  the  one  steady  point  ( x  =  —1)  to  the 
other  (x  =  1).  The  observations  are  also  supported  by  the  comparison  in  Figures  10-11. 

Example  8.3  (Revisted)  Next,  we  also  apply  our  mixed  Gaussian  method  to  Example 
8.3.  We  employ  the  two  Gaussian  distribution  starting  from 

3$  =  [3  x  10s,  2  x  IQ”4,  3  x  10-5f  and  4|o  =  t3  x  10'5>  2  x  10“4^  1(rlF 
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Figure  8:  The  comparison  of  signal,  estimates  of  two  mixed  Gaussian  filter  and  optimal  filter 


Time 


Figure  9:  The  weights  of  two  mixed  Gaussian  filter 
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Figure  10:  The  comparison  of  the  probability  density  functions  of  mixed  Gaussian  filter  and 
optimal  filter  at  t  =  1 


Figure  11:  The  comparison  of  the  probability  density  functions  of  mixed  Gaussian  filter  and 
optimal  filter  at  t  —  2.75 
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with  the  same  covariance 


'  106  0  O' 

P0|o  =  0  4  x  106  0 

0  0  10“4  _ 

respectively  and  the  initial  weights  are  Oq  1  =  Oq2*  =  |.  We  demonstrate  the  performance  of 
the  mixed  Gaussian  filter  with  At  =  0.1  in  Figures  12-15.  We  see  that  the  mixed  Gaussian 
filter  based  on  the  update  formula  (5.5)-(5.6)  captures  the  signal  process  very  well.  That 
is,  the  weights  of  the  two  mixed  Gaussian  filters  shown  in  Figure  15  effectively  pick  up  the 
combined  performance  of  the  Gaussian  filters:  GHF1  and  GHF2.  More  detailed  discussions 
and  numerical  testings  of  the  mixed  Gaussian  filter  will  be  presented  in  the  forthcoming 
paper. 

9  Conclusions 

The  paper  presents  the  systematic  formulation  of  Gaussian  filters  and  mixed  Gaussian  filters 
based  on  the  efficient  numerical  intergation  of  the  Bayesian  formula  for  optimal  recursive 
filter.  Based  on  our  formulation  we  develop  the  two  filter  algorithms,  namely,  the  Gauss- 
Hermite  filter  (GHF)  and  the  central  difference  filter  (CDF).  We  demonstrated  the  feasibility 
of  our  proposed  filter  algorithms  for  testing  nonlinear  filtering  problems.  Our  numerical 
results  indicate  that  both  the  Gauss- Hermite  filter  and  the  central  difference  filter  have 
superior  performance  to  the  filter  of  Julier-Uhlmann  and  the  extended  Kalman  filter.  We 
also  proposed  the  new  update  rules  for  the  Gauss  sum  filters  and  show  that  they  can  perform 
near  optimally. 
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Figure  12:  The  comparison  of  altitude  errors  of  each  Gaussian  filter 
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Figure  13:  The  comparison  of  velocity  errors  of  each  Gaussian  filter 
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